library(tidyverse)
library(qualtRics)
library(cjoint)
library(estimatr)
library(haven)
library(jtools)
library(ggthemes)
library(gridExtra)
library(texreg)
library(ggrepel)
library(ggpubr)
library(patchwork)


# Steps:
# Update figures + text as needed

####################################################
#  Prepare  Data
####################################################

# Load survey responses 
load("data/combined_waves.Rdata")

# Subset to respondents in both waves with consistent age responses
sdata <- subset(sdata,both_waves==1 & abs(w2_age-age)<=1) 

# Load Wave 2 conjoint
cdata <- read.qualtrics("data/wave2_numeric.csv",responses=c("Q15.1","Q160","Q161","Q15.3","Q15.5"),
                        new.format=T,respondentID="ResponseId",covariates=c("PID")) %>% mutate(PID=as.numeric(PID))

cdata <- subset(cdata,!is.na(selected)) %>%  left_join(sdata,by='PID') 
cdata <- cdata %>% select(-Condizione.di.salute.rowpos,-Condizione.economica.rowpos,-Condizione.lavorativa.rowpos,-`Eta.(anni).rowpos`,
                          -Istruzione.rowpos,-Luogo.di.nascita.rowpos,-Partito.politico.rowpos,-Sesso.rowpos)

# Load Wave 1 conjoint: 
cdata_w1 <- read.qualtrics("data/wave1_numeric.csv",responses=c("Q15.1","Q15.3","Q15.5","Q15.7"),new.format=T,respondentID="ResponseId",covariates=c("PID")) %>%
  mutate(PID=as.numeric(PID)) %>% filter(Luogo.di.nascita!="Africa" & !is.na(selected))
cdata_w1 <- cdata_w1 %>%  left_join(sdata,by="PID") %>% filter(!is.na(both_waves))
cdata_w1 <- cdata_w1 %>% select(-Condizione.di.salute.rowpos,-Condizione.economica.rowpos,-`Età.(anni).rowpos`,
                                -Istruzione.rowpos,-Luogo.di.nascita.rowpos,-Partito.politico.rowpos,-Sesso.rowpos) 

# English conversion
source("translation.R")

####################################################
#  Main Text Figures
####################################################

baselines <- list()
baselines$Place.of.Birth <- "Northern Italy"
baselines$Education <- "Lower Secondary"
baselines$Economic.Status <- "Middle class"
baselines$Political.ID <- "Movimento 5 Stelle"
baselines$Age <- "Under 30"
baselines$Occupation <- "Homeworker"
baselines$Gender <- "Woman"

baselines_w1 <- list()
baselines_w1$Place.of.Birth <- "Northern Italy"
baselines_w1$Education <- "Lower Secondary"
baselines_w1$Economic.Status <- "Middle class"
baselines_w1$Political.ID <- "Movimento 5 Stelle"
baselines_w1$Age <- "Under 30"
baselines_w1$Gender <- "Woman"

x <- theme_minimal() + theme(axis.text.y = element_text(hjust = (0)),legend.position = 'none',text=element_text(size=11))

# Figure 1
results <- amce(selected ~ Health.Status + Economic.Status + Education + Occupation + Place.of.Birth + Political.ID + Age + Gender, data=cdata,cluster=TRUE, respondent.id="respondent",baselines=baselines)
pdf(file = "./fig1.pdf",  width = 9, height = 7)
plot(results,main="",plot.theme=x)
dev.off()

# Figure 2
results2 <- amce(selected ~ Health.Status + Economic.Status + Education + Place.of.Birth + Political.ID + Age + Gender, data=cdata_w1,cluster=TRUE, respondent.id="respondent",baselines=baselines_w1)
pdf(file = "./fig2.pdf",  width = 9, height = 7)
plot(results2,main="",plot.theme=x)
dev.off()

# Setup for Figures 3 and 4
source("functions.R")
  
# Figure 3
pdf(file = "./fig3.pdf",  width = 9, height = 4.5)
vis_results("left_right_binary", c("Right Wing", "Left Wing"),wave=2,fig=2)
dev.off()

# Figure 4
pdf(file = "./fig4.pdf",  width = 9, height = 4.5)
vis_results("left_right_binary", c("Right Wing", "Left Wing"),wave=1,fig=2)
dev.off()

# Setup for Figure 5
sdata2 <- sdata
sdata2 <- sdata2 %>% ungroup() %>% rowwise() %>% mutate(
  outgroup_w1 = case_when(
    votechoice_current == 1 ~ mean(c(affect_pd,affect_fdi,affect_fi,affect_lega,affect_ivr,na.rm=T)),
    votechoice_current == 2 ~ mean(c(affect_pd,affect_fdi,affect_fi,affect_5star,affect_ivr,na.rm=T)),
    votechoice_current == 3 ~ mean(c(affect_lega,affect_fdi,affect_pd,affect_5star,affect_ivr,na.rm=T)),
    votechoice_current == 4  ~ mean(c(affect_lega,affect_fdi,affect_fi,affect_5star,affect_ivr,na.rm=T)),
    votechoice_current == 5 ~ mean(c(affect_lega,affect_fi,affect_pd,affect_5star,affect_ivr,na.rm=T)),
    votechoice_current == 6 ~ mean(c(affect_lega,affect_fdi,affect_pd,affect_5star,affect_fi,na.rm=T)),
    votechoice_current == 9 ~ mean(c(affect_lega,affect_fdi,affect_pd,affect_5star,affect_fi,affect_ivr,na.rm=T)),
    votechoice_current == 10 ~ mean(c(affect_lega,affect_fdi,affect_5star,affect_fi,affect_ivr,na.rm=T))
  ),
  outgroup_w2 = case_when(
    votechoice_current == 1 ~ mean(c(w2_affect_pd,w2_affect_fdi,w2_affect_fi,w2_affect_lega,w2_affect_ivr,na.rm=T)),
    votechoice_current == 2 ~ mean(c(w2_affect_pd,w2_affect_fdi,w2_affect_fi,w2_affect_5star,w2_affect_ivr,na.rm=T)),
    votechoice_current == 3 ~ mean(c(w2_affect_lega,w2_affect_fdi,w2_affect_pd,w2_affect_5star,w2_affect_ivr,na.rm=T)),
    votechoice_current == 4 ~ mean(c(w2_affect_lega,w2_affect_fdi,w2_affect_fi,w2_affect_5star,w2_affect_ivr,na.rm=T)),
    votechoice_current == 5 ~ mean(c(w2_affect_lega,w2_affect_fi,w2_affect_pd,w2_affect_5star,w2_affect_ivr,na.rm=T)),
    votechoice_current == 6 ~ mean(c(w2_affect_lega,w2_affect_fdi,w2_affect_pd,w2_affect_5star,w2_affect_fi,na.rm=T)),
    votechoice_current == 9 ~ mean(c(w2_affect_lega,w2_affect_fdi,w2_affect_pd,w2_affect_5star,w2_affect_fi,w2_affect_ivr,na.rm=T)),
    votechoice_current == 10 ~ mean(c(w2_affect_lega,w2_affect_fdi,w2_affect_5star,w2_affect_fi,w2_affect_ivr,na.rm=T))    
  ),  
  diff_outgroup = -1*(outgroup_w2 - outgroup_w1)
)

sdata2 <- sdata2 %>% mutate(
  diff_ui_delta = (w2_ui_italy - ui_italy) - (w2_ui_immigrants - ui_immigrants),
  diff_covid = ifelse(covid_death==1 | covid_serious==1 | covid_ill == 1,1,ifelse(!is.na(covid_death),0,NA))
)

sdata2$diff_outgroup[sdata2$votechoice_current != sdata2$w2_votechoice_current] <- NA

# Figure 5
covars <- "+ gender + factor(age) +edu4 +edu5 + edu6 + edu7 + leftright + 
factor(region) + non_citizen + employed_2019 + factor(sector) + factor(political_interest)"

x <- "diff_outgroup"
out1a <- lm(formula(paste0(x,"~ large_income_instability  +  diff_covid")),data=sdata2)
out1b <- lm(formula(paste0(x,"~ large_income_instability +  diff_covid",covars)),data=sdata2)
out1 <- plot_coefs(out1a,out1b,inner_ci_level=.9,coefs=c("Loss of Income" = "large_income_instability","COVID Exposure" =  "diff_covid"),
                   model.names=c("No Covariates","With Covariates"),legend.title = "",colors=c("darkgray","black")) + 
  ggtitle("Affective Polarization")+ xlab("")   + theme_bw() + ylab("")+ theme(panel.grid.major.y = element_blank())+ 
  theme(panel.grid.minor.x = element_blank())+scale_x_continuous(limits=c(-.7,.7))

x  <- "diff_ui_delta"
out2a <- lm(formula(paste0(x,"~ large_income_instability +  diff_covid")),data=sdata2)
out2b <- lm(formula(paste0(x,"~ large_income_instability  +  diff_covid",covars)),data=sdata2)
out2 <- plot_coefs(out2a,out2b,inner_ci_level=.9,
                   coefs=c("Loss of Income" = "large_income_instability","COVID Exposure" =  "diff_covid"),
                   model.names=c("No Covariates","With Covariates"),legend.title = "",colors=c("darkgray","black")) + 
  ggtitle("Welfare Chauvinism") + theme_bw() + theme(axis.text.y=element_blank()) + xlab("") + ylab("") + 
  theme(panel.grid.major.y = element_blank())+ theme(panel.grid.minor.x = element_blank()) +
  scale_x_continuous(limits=c(-.7,.7))

pdf(file = "./fig5.pdf",  width = 9, height = 4.5)
out1 + out2 + plot_layout(guides = 'collect') & theme(legend.position='bottom')
dev.off()

# Table A3
texreg(list(
  out1a,out1b,out2a,out2b), 
  caption.above=TRUE,
  include.rsquared = FALSE,
  include.adjrs = FALSE,
  include.groups = FALSE,
  stars=numeric(0),
  custom.coef.names=c("Intercept","Income Shock","Health Shock","Female","Age: 18-24","Age: 25-34","Age: 35-44","Age: 45-54","Age: 55-64","Age: 65+",
                      "Education: Vocational","Education: Preparatory","Education: College","Education: Advanced",
                      "Left-Right Ideology","Region 2","Region 3","Region 4","Region 5",
                      "Region 6","Region 7", "Non-citizen","Employed in 2019","Sector - Manufacturing","Sector - Utilities",
                      "Sector - Construction","Sector - Retail","Sector - Finance","Sector - Real Estate","Sector - Hospitality","Sector - Transport","Sector - IT","Sector - Education","Sector - Health",
                      "Sector - Public Admin","Sector - Other","Political Interest - Interested","Political Interest - Little interest","Political Interest - Not at all Interested"),
  digits=3, scalebox=0.75, booktabs=TRUE,  use.packages = FALSE,float.pos="!h")



####################################################
#  Appendix Items
####################################################

#### Appendix A

# Load survey responses 
load("data/combined_waves.Rdata")

# Table A1
sdata$both_waves[sdata$both_waves == 0] <- "W1"
sdata$both_waves[sdata$both_waves == 1] <- "W2"

balance <- sdata %>%
  select(both_waves, female, age1, age2, age3, age4, age5,age6, edu1, edu4,edu5,edu6,edu7) %>%
  gather(key=variable, value=value, -both_waves) %>%
  group_by(both_waves,variable) %>%
  summarise(value=list(value)) %>%
  spread(both_waves,value) %>%
  group_by(variable) %>%
  mutate(mean_w1=mean(unlist(W1),na.rm=T),mean_w2=mean(unlist(W2),na.rm=T),t_value=-1*t.test(unlist(W1),unlist(W2))$statistic)

balance <- balance %>% ungroup() %>% select(-W1,-W2,-variable)

balance <- as.data.frame(balance)
row.names(balance) <- c("Age: 18-24","Age: 25-34","Age: 35-44","Age: 45-54","Age: 55-64","Age: 65+",
                        "Education: Secondary","Education: Vocational","Education: Preparatory","Education: College","Education: Advanced","Female")

knitr::kable(balance,digits=3,col.names=c("Wave 1 Mean","Wave 2 Mean","t-statistic"),row.names=T)




# Create new binary variables

cdata <- cdata %>% mutate(willvote_pd=ifelse(w2_votechoice_current==4,1,0))

# Figure A1
pdf(file = "./appendix_figa1.pdf",  width = 9, height = 8)
vis_results("left_right_binary", c("Right Wing", "Left Wing"),wave=2,fig=1)
dev.off()

# Figure A2
pdf(file = "./appendix_figa2_a.pdf",  width = 9, height = 8)
vis_results("willvote_pd", c("Vote Intention != PD", "Vote Intention PD"),wave=2,fig=1)
dev.off()
pdf(file = "./appendix_figa2_b.pdf",  width = 9, height = 4.5)
vis_results("willvote_pd", c("Vote Intention != PD", "Vote Intention PD"),wave=2,fig=2)
dev.off()

# Figure A3
pdf(file = "./appendix_figa3.pdf",  width = 9, height = 8)
vis_results("left_right_binary", c("Right Wing", "Left Wing"),wave=1,fig=1)
dev.off()

# Table A2

with(sdata2,prop.table(table(large_income_instability,diff_covid)))

# Table A3

  # See main text section

# Figure A4
pdf(file = "./appendix_figa4.pdf",  width = 12, height = 7)
vis_results("large_income_instability", c("No Significant Income Loss", "Significant Income Loss"),wave=2,fig=2)
dev.off()

#### Appendix B


# Remove middle tercile for GAL/TAN and Econ Left/Right
cdata$galtan_q3[cdata$galtan_q3 == 1] <-NA
cdata$galtan_q3[cdata$galtan_q3 == 2] <-1
cdata$lrecon_q3[cdata$lrecon_q3 == 1] <-NA
cdata$lrecon_q3[cdata$lrecon_q3 == 2] <-1

# Figure B1
pdf(file = "./appendix_figb1_a.pdf",  width = 9, height = 4.5)
vis_results("galtan_q3", c("GAL", "TAN"),wave=2,fig=2)
dev.off()

pdf(file = "./appendix_figb1_b.pdf",  width = 9, height = 4.5)
vis_results("lrecon_q3", c("Left economic","Right Economic"),wave=2,fig=2)
dev.off()

# Figure B2
pdf(file = "./appendix_figb2.pdf",  width = 9, height = 8)
vis_results("galtan_q3", c("GAL", "TAN"),wave=2,fig=1)
dev.off()

# Figure B3
pdf(file = "./appendix_figb3.pdf",  width = 9, height = 8)
vis_results("lrecon_q3", c("Left economic","Right Economic"),wave=2,fig=1)
dev.off()

#### Appendix C

# Load survey responses 
load("data/combined_waves.Rdata")

# Figure C1
dat <- sdata %>% mutate(left_right_binary2=ifelse(left_wing==1,"Left",ifelse(right_wing==1,"Right",NA))) %>%
  select(covid,left_right_binary2) %>% na.omit() %>%
  count(left_right_binary2, covid) %>%
  group_by(left_right_binary2) %>%
  mutate(percent = n / sum(n),
         error = sqrt((percent * (1-percent))/n)) %>% arrange(desc(percent))

pdf(file = "./appendix_figc1.pdf",  width = 9, height = 4.5)
ggplot(dat, aes(x=reorder(covid,-percent), y = percent, fill = as.factor(left_right_binary2)))+ 
  geom_bar(position = "dodge", stat = "identity") + 
  geom_errorbar(stat="identity",aes(ymin=percent-error,ymax=percent+error),position=position_dodge(0.9, preserve = "single"),width=.15) + 
  guides(fill=guide_legend("Self-identified Ideology")) + xlab("Covid Exposure") + ylab("Proportion") + theme_bw()
dev.off()

# Table C1
round(prop.table(table(sdata$w2_vaccine_intentions)),3)

# Figure C2
pdf(file = "./appendix_figc2.pdf",  width = 9, height = 7)
ggplot(sdata,aes(x=leftright,y=vaccine_recode)) + geom_smooth() + ggtitle("Intention to take Vaccine") +
  ylab("Percentage of Respondents") + xlab("Left (0) to Right (10) Self-Placement") 
dev.off()

# Figure C3
pdf(file = "./appendix_figc3.pdf",  width = 9, height = 7)
cdata2 <- cdata %>% filter(vaccine_recode==1)
results <- amce(selected ~ Health.Status + Economic.Status + Education + Occupation + Place.of.Birth + Political.ID + Age + Gender, 
                data=cdata2,cluster=TRUE, respondent.id="respondent",baselines=baselines)
plot(results,main="",plot.theme=x)
dev.off()
